
#################################################################################
# File Name:  	main.R 
# Paper:			  Beliefs about Climate Action Consequences under Weak Global Institutions: #               Sectors, Home Bias, and International Embeddedness
# Author: 		  Patrick Bayer (Strathclyde) and Federica Genovese (Essex)
#
# Purpose:      Main results
# Data input:   "./gep_data.csv"
# Date:         26 September 2020
#
#################################################################################

rm(list=ls())

# Load packages
library(ggplot2)
library(dplyr)
library(ggpubr)
library(reshape2)


# Load data
f<-read.csv("gep_data.csv")


# Drop incomplete responses
f <- f[f$Progress==100,]

# Exclude those who failed the attention check
f <- f[f$attentioncheck==3 & is.na(f$attentioncheck)==FALSE,]

# Exclude pacers who were faster than 3 minutes 
f<- f[f$Durationinseconds>180,]
mean(f$Durationinseconds)/60

# Clean data (Don't know responses)
f$devep_ccperformance[f$devep_ccperformance==5] <- NA
f$deving_ccperformance[f$deving_ccperformance==5] <- NA
f$Commonwealth[f$Commonwealth==4] <- NA


# ================================================================
# Descriptives
# ================================================================

# Basic demographics
table(f$gender,useNA="ifany") #51.5% of our sample is female
table(f$age, useNA="ifany") #42% of our sample is 40 or below
table(f$education, useNA="ifany") # 53% has a university degree
summary(f$income) # Median income category of 30-40,000

# Basic political variables
table(f$political_id, useNA="ifany") #46% of our sample is center-left to left
table(f$X2017_election, useNA="ifany") #82% voted 
table(f$Brexit, useNA="ifany") # 60% Remainers, 32% Leavers, 8% Undecided

# Climate change variables
table(f$cc_seriousproblem, useNA="ifany") # 93% thinks climate change is a serious problem
table(f$cc_concern, useNA="ifany") # 79% are somewhat or more concerned about climate change
table(f$heardparisaccord, useNA="ifany") #82% have heard of Paris Agreement, but 60% of those know little about it

# Sectoral variables
table(f$sectorperformance, useNA="ifany") #51% think energy and transport will do poorly
table(f$transportsectorfeel, useNA="ifany") #38% use public transportation
table(f$energysectorfeel, useNA="ifany") #52% understand electricity bill fairly well

# ================================================================
# Variable creation
# ================================================================

# Create indicators for experimental conditions
f$treat_vign[is.na(f$t_control_Page.Submit)==FALSE] <- 0
f$treat_vign[is.na(f$t_t1_domwin_Page.Submi)==FALSE] <- 1
f$treat_vign[is.na(f$t_t2_domlos_Page.Submit)==FALSE] <- 2
f$treat_vign[is.na(f$t_t3_forwin_Page.Submit)==FALSE] <- 3
f$treat_vign[is.na(f$t_t4_forlos_Page.Submit)==FALSE] <- 4

f$treat_vign <- factor(f$treat_vign, levels=c(0,1,2,3,4), labels=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"))

table(f$treat_vign, useNA="ifany")


# ================================================================
# BELIEFS
# ================================================================

# UK winner (strongly loser, slightly loser, slightly winner, slightly winner) 
table(f$vign_beliefUKwins, useNA="ifany")
table(f$vign_beliefUKwins, useNA="ifany")

# Germany winner 
table(f$vign_beliefGerwins, useNA="ifany")
f$vign_beliefGerwins
table(f$vign_beliefGerwins, useNA="ifany")

f$dummy_beliefGerwins <- ifelse(is.na(f$vign_beliefGerwins)==TRUE,NA,ifelse(f$vign_beliefGerwins>2,1,0))
table(f$dummy_beliefGerwins)


# ================================================================
# MAIN ANALYSIS
# The Effect of Winning/Losing Sector Information and Home Bias
# ================================================================

# Belief that UK wins
h1 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: UK wins
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)

# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
pvallist[i-1] <- t.test(vign_beliefUKwins ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")



# Belief that Germany wins
h2 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefGerwins, na.rm = TRUE),
    sd = sd(vign_beliefGerwins, na.rm = TRUE)
  )

# ANOVA: Germany wins
res.aov <- aov(vign_beliefGerwins ~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)


# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_beliefGerwins ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")



# ================================================================
# FIGURE 1: Belief UK/Germany win 
# ================================================================

# Restructure data for better plotting options
df.plot <- cbind(f[,c("treat_vign","vign_beliefUKwins","vign_beliefGerwins")])
colnames(df.plot)[2:3] <- c("Belief: UK wins", "Belief: Germany wins")
g <- melt(df.plot,value.name="outcome")                    


pdf(file="./belief_UKGER_win.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c("Belief: UK wins", "Belief: Germany wins"), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(g, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
            ) +
            facet_grid(~variable, scales="free") +
            scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
            geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
            theme_bw() +
            theme(text = element_text(size=11),
                  axis.text.x = element_text(size=9),
                  legend.position="none"
                  )
dev.off()
            



# ================================================================
# BREXIT SUBGROUP ANALYSIS
# The Effect of International Embeddedness
# ================================================================

f$BrexitBin <- ifelse(f$Brexit==1,1,ifelse(f$Brexit==2,2,NA))
subdf <- f[is.na(f$BrexitBin)==FALSE,]
subdf$BrexitBin <- factor(subdf$BrexitBin, levels=c(1,2), labels=c("Leavers", "Remainers"))

# Beliefs for Leavers
h1 <- group_by(subdf[subdf$BrexitBin=="Leavers",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Leavers
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$BrexitBin=="Leavers",])
summary(res.aov)
TukeyHSD(res.aov)


# p-values for t-test of treatment effects (Leavers)
pvallist <- NA
temp <- subdf[subdf$BrexitBin=="Leavers",]

for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_beliefUKwins ~ treat_vign, data=temp[as.numeric(temp$treat_vign)==1 | as.numeric(temp$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# Beliefs for Remainers
h2 <- group_by(subdf[subdf$BrexitBin=="Remainers",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# ANOVA: Remainers
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$BrexitBin=="Remainers",])
summary(res.aov)
TukeyHSD(res.aov)

# p-values for t-test of treatment effects (Remainers)
pvallist <- NA
temp <- subdf[subdf$BrexitBin=="Remainers",]

for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_beliefUKwins ~ treat_vign, data=temp[as.numeric(temp$treat_vign)==1 | as.numeric(temp$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")



# Compare Leavers vs Remainers
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ BrexitBin, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# FIGURE 2: Brexit subgroup
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","BrexitBin","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")
 

pdf(file="./belief_UK_win_Brexit.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
            ) +
            facet_grid(~variable, scales="free") +
            scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
            geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
            theme_bw() +
            theme(text = element_text(size=11),
              axis.text.x = element_text(size=9),
              legend.position="none"
              )
dev.off()


# ================================================================
#                           END OF FILE
# ================================================================
